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ABSTRACT 

We study the two dimensional Hubbard model by use of the 
ground state algorithm in the Monte Carlo simulation. We employ 
complex wave functions as trial function in order to have a close look 
at properties such as chiral spin order (xSO) and flux phase. For 
half filling, a particle-hole transformation leads to sum rules with 
respect to the Green's functions for a certain choice of a set of wave 
functions. It is then analytically shown that the sum rules lead to the 
absence of the x^O. Upon doping, we are confronted with the sign 
problem, which in our case appears as a "phase problem" due to the 
phase of the Monte Carlo weights. The average of the phase shows 
an exponential decay as a function of inverse temperature similarly 
to that of sign by Loh Jr. et. al. . We compare the numerical results 
with those of exact numerical calculations. 
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§1. Introduction 

The two dimensional Hubbard model has recently attracted much attention 
in view of possibly providing insight into the high-Tg Cu-0 superconductors. In 
the strong U (Coulomb repulsive force) limit at half filling, the model is equiv- 
alent to the antiferromagnetic Heisenberg model, and in good agreement with 
the antiferromagnetic behavior of superconducting materials. When the system is 
doped, however, the ground state properties appear less clear. Concerning ground 
states of the strongly correlated electronic systems, various hypothetical states 

1) 2 3) 4) 

such as spin-liquid , resonating valence bond( RVB) , ' and flux phase, etc. 
are proposed. Of particular interest is a possibility of the violation of P (parity) 
and T (time) invariance in such systems."'^'' From the field theoretical point of 
view also, the relevance of spontaneous breakdown of such symmetries have been 
pointed out, and it is conjectured that a Chern-Simons term in 2+1 dimensional 

5 — 7) 

QED, which is regarded as an effective theory describing the RVB picture, may 
violate the confinement of the fundamental charge ( the charge-spin separation). 
In spite of such intriguing arguments, however, any clear evidence has not so far 
been found. 

The above ideas are proposed based basically upon the mean field approach, 
and therefore the stability of its vacuum against full quantum fluctuations is not 
clear at all. In such situations numerical approach is a suitable method to study 
the strongly correlated electron systems. Among various approaches, Monte Carlo 
simulation is one of the most promising methods. 
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In the present paper we present a Monte Carlo study ' of the two dimen- 
sional Hubbard model/ ""^^ So far similar works have been done. Although our 
approach is basically the same, and the size of lattice is also quite limited, our 
motivation of the study is two fold. One is to use complex wave functions as 
a trial function in the ground state algorithm in order to have a close look at 
the properties stated above, namely, ones of the chiral spin order (xSO) and flux 

4 12) 

phase, ' where the complexity of the expectation values is involved. At half- 
filling, a particle-hole symmetry makes the theory transparent to look at. By use 
of such symmetry, for some choice of a set of wave functions (we call it "dual 
choice"), it is shown that there hold sum rules with respect to the Green's func- 
tion with spin up and that with down. They are of different form for even and odd 
lattice spacings. We analytically show that these sum rules lead to the absence of 

13) 

the expectation value of the xSO parameter at half-filling. 

Another motivation is to study an effect of such wave functions to the so called 
sign problem, which frequently arises in the quantum Monte Carlo simulations. 
A widely used algorithm employs the Hubbard-Stratonovich transformation,^'^^^ 
rewriting the system to the 2-1-1 dimensional classical system of the Ising dynam- 
ical variable. The conventional importance sampling with configurations of this 
variable encounters the negative weight which is originated from the fermionic de- 
terminant. This prevents probabilistic algorithm from being applied to the system. 
A way to circumvent the problem is to use the positive weight and to estimate 
expectation value of operator combined with the average value of the sign. This 
however causes large statistical errors for some cases. 
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The origin of the sign problem in the framework of trial function approach 
is not yet clear, and thus a way to overcome it has not yet been established. 

14) 

Some improvements have been tried by employing optimized wave function. 
The results depend very much on the choice of the wave functions of the trial 
state. It is therefore not trivial at all what to obtain if we employ the complex 
wave function. In our case the sign of the weights are replaced by their phase 
("phase problem") . We paid a special attention to the average of the phase. In 
the present paper, we study, as doped cases, 14 fermion and 10 fermion systems 
on 4x4 lattice. The latter case is almost free from the phase problem, but the 
former receives rather serious effect. A measurement of the x^O for the former 
case unfortunately gets so large noises to extrapolate to the low temperature limit. 
However other physical quantities such as one particle density operator and spin- 
spin correlations are rather in good agreement with the result obtained by exact 
numerical calculations.^^ ^^"^ 

This paper is organized as follows. In the following section, we present the 
formalism. In section 3, the half filled case is analytically discussed based upon a 
particle-hole transformation. Some choice of a set of wave functions is shown to 
lead to a sum rule. The absence of the expectation value of xSO is derived at half 
filling. In section 4, results of numerical simulations are presented for half filled 
and doped cases. Our conclusions and discussion are presented in section 5 . 
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§2 Hubbard model and Monte Carlo algorithm 

In this section, we shortly present the formahsm for fixing notations as well as 
for making the paper self-contained. It is basically based on that of Hirsch^'^^ 
and Imada-Hatsugai.^°^ 

The Hubbard hamiltonian is defined as 

H = -tJ2 4aCj> + h.c. + uJ2 ^n^a (2-1) 

{ij)a i 

where Cjo.(c|^) is an annihilation (creation) operator of an electron with spin a 
(cr = up or down) at site i, and njo- is the number operator at site i of spin 
a. The first term is the hopping term of electron over nearest neighbor sites 
(ij), and the second stands for the on-site repulsive Coulomb potential with the 
strength U{> 0). An algorithm to perform Monte Carlo simulations here is so 
called ground state algorithm for fixed number of fermions. Expectation value of 
a physical observable O is calculated by 

_ ($|Oexp(-/3g)|$) 
^ ^ " p(/3; $) ^ 

where 

p(/3;$) = ($|exp(-/?iy)|$) (2.3) 

is given in terms of an appropriate trial state vector |$). As /? goes to infinity, 
(2.2) gives an expectation value for true ground state IV'o), unless the trial function 
is orthogonal to the true ground state; 

exp(-/3if)|$) ^ exp(-/3i?o)|Vo)(^o|$) (2.4) 
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For evaluating (2.2) we make Monte Carlo simulations. Introducing discrete pseu- 
dotime with interval Ar, p(/?; is rewritten as 

= ($|(exp(-ATiy))^|$) (2.5) 

where At — (3/L and L denotes the Trotter size, i.e., the number of time slices. 
Separating the hopping term Hq and the self-interaction term Hi of the hamilto- 

21 22) 

nian, one introduces Trotter-Suzuki formula. 



exp(-ATi?) = exp(-^iyo) exp(-ATi?i) exp(-^iyo) + ^((At)^) (2.6) 



To deal with the quartic interaction term Hi of the fermion field, one introduces 
the Hubbard-Stratonovich (H-S) transformation of the Ising type for each site. 



exp(-cn^n4) = ^ ^ exp (2as(nt - n;) - |(nt + n;)) (2.7) 

s=±l 

where a = tanh""*^ •\/tanh(c/4) for a constant c. Eqs.(2.5), (2.6) and (2.7) lead to 

L 

p(A$)= ^($1 X[{woWi{{su})wo)m (2.8) 

{Sil} 1 = 1 

where the summation is taken over the H-S variable su on the 2+1 dimensional 
lattice. Quantities wq and wi{{sii}) are defined as 



wq =exp(— Ari^o/S) 

Art/, (2.9) 



w^i({sii}) =n(^ cxp[2a[7Sii (nt - n|) - ^^(nt + ^4)]) 



where ajj = tanh ^ iytanh(Art/'/4), and A?" is the number of lattice sites. 
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We assume a factorization of spin up and down sectors for the trial state, 

|$) = |$^)|$;) (2.10) 
and take wave functions of single particle basis for Mf fermions for each spin sector 

Mf N 

i*'^)=n(E('^-)-^4)io) (2.11) 

m i 

where {(j)a)mi{i — 1, N;m = 1, Mf) is a Mf x N matrix. In this paper we 
allow {(f)a)mi to be complex in general. Using (2.11), p(/?;$) (2.8), is explicitly- 
represented in a matrix form 

p(/3;$) = J] W^(/3;s) 

{su} (2.12) 
WiP;s) =W^iP;s)WiiP;s) 

where Wf^{P; s) is given by a determinant of a product of matrices 

L 

W,{l3;s) = det[(l>l J](Mo Mi,(0 Mo) ] (2.13) 

Matrix Mq is given by 

Mo =exp(-i^) 

_ ( —ATt/2 nearest neighbor (i,j) 
'■^ 1 otherwise 

while Ml, which depends upon the H-S variable, is a diagonal matrix with element 

{Mi^{l))u = exp[aauSii - At/2] (2.15) 
for cr = +1(— 1) for spin up (down). 

8 



(2.14) 



According to Imada-Hatsugai, single particle Green's function at the time 



slice I 

{La 

(njn v.. = — 

(L,(0|i?.(0) 



( ^-(0 h = \ r mio m (2-16) 



is expressed as 

G^{1) = R^{l)g^{l)Ll{l) (2.17) 

where 

^,(0 = (4 (Oi?,(0)-' (2.18) 

Matrix Ra{l) is N x Mf matrix, and Ll{l) is Mf x N one; 

L 

Rail) =M^a{l)Mo \[ {MoM^,{l)Mo)^, 



1=1+1 

l-l 



(2.19) 



=4 n ( MoM^,il)Mo ) Mo 
1=1 

Finally, expectation value of operator with spin a for the ground state is 

averaged at each time slice, and given by 

, _ r E{.} det ( L.(Ot R^(l) ) 

^1 ^ det {LaiiyOg Rail)) 1 ^^-^^^ 

We employ the heat-bath method for updating the H-S variables. W, however, 
may not be real in some cases, though p(/3; $) is positive semi-definite. These 
cases lead one to so called sign problem. In our case it is replaced by that of the 
phase, which prevents us from doing importance sampling. Conventional way out 
of it is to define a positive measure 

\W\ 



+ 



Es\w\ 
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(2.21) 



and evaluate (2.20) in such a way as 



^ l^ det(L.(Oto^i?.(0) 1 



^^1 ^ det(L.(/)^Q.i?.(/)) ^,,^„^^^,^^,,,„^^ 

(2.22) 

where = |VF|e*^'", and (•)+ is the expectation value of • with the positive 
measure P+. 

§3 Dual choice of wave function and sum rule 

The hamiltonian of the system respects a symmetry under particle-hole trans- 
formation 

cl^fj%-^, (3.1) 

where di-^r denotes the annihilation operator of a hole at site i with spin —a . The 
sign factor rj'^ = (—1)* is positive (negative) for even (odd) sites. In this section 
we will show the Green's functions satisfy sum rules when the system is half-filled, 
and "dual choice" of the wave function is made. The sum rule is satisfied when 
the particle-hole symmetry holds. The sum rules are 

f _ / Sij - 9*, for \i - j\ =even 



for|z-j|=odd 

where and gij denote the spin up and down Green's function in the configura- 
tion space, respectively, and * is the complex conjugation. 

S-1 particle-hole transformation 
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Particle-hole (p-h) transformation gives the correspondence 

|0) < ^ |2) \la) < ^ \la) |2) < ^ |0) 

at each site, where |0), \ la) and |2)denotes zero, single particle with spin a and 
two particles (spin up and down), while |0), |icr) and |2)denotes zero hole, single 
hole with spin a and two holes, respectively. Note that hole states are denoted 
with tilde. 

The creation operator of a particle with quantum number m (e.g., momentum) 
is given as 

j 

where (prnj is the wave function. Suffices m and j are for example abbreviated as 
m = k = {ki,k2),j = n = (ni,n2) for the momentum {ki,k2) and the position 
(ni,n2). Let us see how states transform under the p-h transformation. For 
illustration we consider one-dimensional half-filled {Mf = 2) system with AT = 4 
(4 site system). Spin a state is 

where mi and m2 are quantum numbers of the two particles with spin a. Under 
the p-h transformation, this state goes to 

X]0miii</>m2i2^?ii^i2t^ii-aC^i2-a4-a4-<74-<74-a4a4<74a4<7|O) (3-5) 

We pay an attention to spin —a quantities, and then rewrite 

dUA-a4-A-a = E ^hhlslA-A-A-A-a (3-6) 

11121314 
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where Zi, .., I4 denote the sites and Si^i^i^i^ stands for the Levi-Civita symbol. The 
quantities with —a in (3.5) is 

jij2lihl3h 

= Yl ^rmh(l>m2l2VliVl2£lil2l3Udl-adu-a\^) (3.7) 

11121314 

— ~ 2^ ^ V ^m3j34'm4j4''lj3Vj4^mim2m3m4dj^-crdj4-cr\^) 
j3j4m3m4 

where A = det0 and (f)mj = 4^mjVj- Here we assumed the orthonormal condition 
of wave functions 

j j 
Y ^^j^mj' = Y ^^J^*mj'VjVj = Sjj' 



(3.8) 



and used an identity 

^iii2j3i4'?^miii0m2i2'?^rn3j3 0m4j4 ~ £mim2m3m4^ (3-9) 

From Eq.(3.8) and (3.9) it foUows that 

^jlj2j3j4 0m3i3 07714^4 ~ ^mijl^m2j2^^i^^'>^3m4 

A (3.10) 

which leads to (3.7). Note that the quantum numbers ms, m4 on the r.h.s. of (3.7) 
turn out to be different from the ones (fixed) mi, m2 on the l.h.s. of (3.5). It is due 
to the factor £rnim2m3m4) or the Pauli's principle. The consequence of (3.9) is ; this 
state is similar to the original one (l.h.s. of (3.4)) if we replace (i) c''' by d\ (ii) a by 
—a, (iii) (f)mj by ^j^^- = 4>^jr]j , where rh is the coset of m, namely a combination 
which is completely different from the original set of quantum number m. Having 
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in mind the fact that the hamiltonian is symmetric under p-h transformation, the 
original Green's function of spin o based on the wave function is transformed 
to that of spin — cr based on 4>iny 
3-2 dual choice 

Consider a set of wave function {4>mj} which defines (3.6) of single particle 
basis. We learned in the previous sub-section that the corresponding hole basis is 
i^fhjVj}- When the set {0mj^j} coincides with the set {4>mj}, that is, when the 
hole basis wave function in (3.9) is simply complex conjugate of the particle basis 
wave function, we call the choice of the wave functions {(f>mj} " dual choice ". Let 
us take the previous example of the one dimensional 4-site case. We can take (prnj 
as 

= 1 (j)2j = Cexp(i27rn/4) (3.11) 

where n is an integer standing for the location j, and C is a normalization constant. 
Then can be taken as 

(plj = Cexp(-i27rn/4) = Cexp(i47rn/4). (3.12) 

We see that 

(PljVj = <l>*2j (f>ljVj = (Plj, (3.13) 
and it shows the choice (3.11) is a "dual choice". 
3-3 sum rule 

Since the hamiltonian is invariant under the p-h transformation, we drop for 
simplicity the Boltzmann factor exp(— /Ji?) in the following arguments. Due to 
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the dual choice of wave function, matrix element of c's is transformed to that of 
d's with the same set of wave functions except for complex conjugation. 

^ , , ~ (3.14) 

jlj2j[j2 

The Green's function is changed to 

Gaiji^ffy ~{^\'^rn2aCmiaCiac}jcr^mia^m2(j\^) 

=Sij - r]ir]jG-crji{4>*) 

(3.15) 

where the reality of G is used. We therefore have 



h={^'^~'^' (3.16) 

' 5f*j for \ =odd ^ 



where /y and denote the spin up and down Green's function in the configura- 
tion space. 

3-4 general case 

We now turn to general case where the lattice size N is Ni x N2 {Ni and 
A^2 stand for the number of lattice sites in the first and second direction) and the 
fermion number of each spin M = N/2 (N =even). M fermion state is 

^rni,a'^rn2a- ' ' ' '4nMcr\^) ~ 'i^rmji ' ' ' 'PmM jM^j'i_<j^j2a ' ' ' ^\Mcr\^) ' (3.17) 
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It changes to 



3i32---3m 

1 

iV! 



(3.18) 



Concerning the spin —a quantities, we have 

nCmMI 



^1 0miii ■ ■ ■ ^mMiMViiVi2 ' ' ' V 

ii ...ijv 



^M 



'^^iii2---iMiM+i---iN(^iM+i ' ' ' ^In^^^ 
{j}mM+i---mN 



(3.19) 



where A = det^, (f)mj = (l>mjVjj mi...mM are fixed. For dual choice, we have 

^mM+i3M+i ' ' ' ^mN3N ~ ^mijM+i ' ' ' ^ruMjN (3.20) 

So (3.20) teUs us that the transformed wave functions can be given by the original 
wave functions. Only the difference from the original wave function is that they 
appear in complex conjugate form. We assume |A| = 1 because of the unitarity 
of the wave function. Along with the similar argument given in §3-3, we reach 

f = i ~ ^""^^ spacing 

1 3ji spacing ^ ' ' 

3-5 x^O vs. sum rule 

In order to calculate the chiral spin order parameter (xSO),^^ six point func- 
tion is necessary. Six point function with the same spin operators is decomposed 
into triple products of two point function (Green's function); 

(C1C1C2C2C3C3) —G11G22GSS ~ G12G21GS3 + G12G31G23 

(3.22) 

— G11G23G32 + G1SG21G32 — G13G22G31 
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The xSO is given by 



-^123 = {^1 ■ {^2X^3)) 



(3.23) 



where (Tj is Pauh matrices sitting at site i, and 



PI123 = 



(Xl2X23X3l) = (cLc2a4a'C3a'4^„Ci^") 



(3.24) 



As is shown in the appendix, (P/123 — -P/132) is purely real provided that the sum 
rule (3.21) holds. So a consequence of the sum rule is that the real part of the 



becomes zero. 

This holds only for the half-filling case and does not tell anything about doped 
cases. In the following section we will explicitly see the above properties by making 
numerical simulations. 

§4. Monte Carlo simulations 

In the present paper we take 4x4 spatial lattice, and the Trotter size L is taken 
to be 100 through 200, so that At is sufficiently small ( At < 0.1). The hopping 
parameter t is kept to be 1, and U is taken to be 4.0. The wave function {0rnO is on 
the free particle basis in the complex form exp(i^k ■ n). We choose particle states 
with momentum k so as to fill the energy states in order from the lowest to higher 
ones. We make a dual choice of section 3. For half-filling, the lower-half energy 
states are occupied by {(f)mi} and the dual counterparts {ipmi} correspond to the 



xso 



ReEi23 = 2Im(PZi23 - ^^^132) 



(3.25) 
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upper-half. Setting {cpmi} = exp[i27rk • n/4], our choice is k = (^1,^2) = (0,0) 
for (f>m = (l>i, (1, 0) for cf>2, (-1, 0) for 03, (0, 1) for (f>^, (0, -1) for 05, (1, 1) for 06 
and (1, —1) for 07. As for 08, we choose a hnear combination of the two states 
(2, 0) and (0, 2). Dual counterparts are chosen in such a way that {ki,k2) = (2, 2) 
for V'l, (-1, 2) for ^^2, (1, 2) for V'a, (2, -1) for ^^4, (2, 1) for V's, (-1, -1) for V'e, 
(—1,1) for 7/^7 and i/jg = 08- In Fig.l, we show them in the momentum space. 
For spin up and down sectors we take the same wave functions {0mi}, and m 
runs from 1 to Mf. We performed from 3000 to 18000 sweeps depending on M/ 
(the number of fermions with each spin a (2.11)) and p. The statistical errors are 
estimated by use of the block method with each block size being 500 to 1000 after 
1000 warming-up sweeps. 

4.1 half filled case (Mf ^8) 

For the half filled case, no sign problem is involved, since the average of phase 
is purely real and is unity for any configurations occurred in the algorithm. We 
define symmetrized two point functions 

1 1 

fm = + fji) 9{ii} = ^{9ij + 9ji) (4.1) 

Spin average of the symmetrized functions satisfy 

^^U{i3} + 9{ij}) = for |i - j I = even 7^ 

^e{f{u}+9{ii})^l-Q (4.2) 

I"^(/{ii} + 9{ij}) = for |i - j| = odd 
In Fig. 2a, we show results of the Green's function, and compare with those of 

the Lanczos's method.^^'^^^ The average values of the real part for even spacings 
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(K ~i| 7^ 0) exactly zero in agreement with the argument in section 3. Fig.2b 
is spin-spin correlations. It shows clear antiferromagnetic correlations and agrees 
very well with the Lanczos's one. The numerical figures for the Green's functions 
and the spin-spin correlations are listed in Table I. 

The real part of xSO turns out exactly vanishing as discussed in the previous 
section. We looked also at the imaginary part, and it is consistent with zero within 
error. 

4-2 sign problem 

For simulations of doped cases one is confronted with the sign problem. In our 
case, for the complex trial wave function, W is generally complex, although p(/3; 
(2.3) is real. As stated in section 2, one generates a sequence of configurations 

23) 

with a probability, 

P.W = (4.3) 

P+ 

where 

p+ = j2ms)\ (4.4) 

s 

The average of the phase of W, which appears in (2.22), is then given by 

< e''- >+= p/p+. (4.5) 
For large P, p is dominated by the ground state Ii/jq) with energy Eq 

p~exp(-/?£;o)| < Vo|*> 1^ (4.6), 
18 



while p+ looks like 

p+ ~ exp(-/?E+)| < V+ol* > I', (4.7) 

where E+ and ip+o are analogue of Eq and ipo for p+, respectively. Therefore the 
average of the phase is real, and behaves in a similar manner to that of the sign 
in the real trial wave function case; 

(e^^» )+ = p/p+ ~ exp{-(3AE) (4.8) 

where AE = Eq — E+. For AE = 0, the average of sign converges to some finite 
value. In this case, one has no difficulty with the phase, and one may neglect the 
phase for evaluating the average of physical quantities. For AE' 7^ 0, on the other 
hand, we encounter so called the sign problem, that is, {e^^^) becomes vanishing 
in our case. It apparently gives meaningless result for the average of any physical 
operator (2.22) and anticipates, at the same time, large errors. If, however, the 
numerator in (2.22) behaves in a exponential manner similarly to the denominator, 
it is still expected that the average may converge to some finite value, and its 
errors could be tamable. This depends upon the fillingness and operators to be 
calculated. In the following subsections we show the results concerned with this 
issue for M/=7 and M/=5 (closed shell). 

4-3 Hole doped cases Mf — 7 

The average of the real part of the phase is shown in Fig. 3 and Table II. In 
the region we have calculated {P = 2.0 to 5.0), an exponential behavior is clearly 
observed. The AE read from the slope is 1.01(5). This value is somewhat larger 
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than the one obtained by Loh et al. The average < O with the measure 
(2.21) for Green's functions (O = Gij) and for spin correlations are determined 
in a high precision, while the real part of the O — xSO parameter is somewhat 
noisy and looks unstable as (3 varies. This difference of the fluctuations in < O >+ 
results in the difference of weighted averages (2.22) among the Green's function, 
spin-spin correlations and the xSO. 

According to the formula (2.22), the average of observables is proportional 
to the product of observables and the phase in the numerator. We have found 
that due to the steadiness of behaviors of the Green's function Gij and spin-spin 
correlations with respect to a variation of /?, their product shows an exponential 
fall-off with almost the same slope as the denominator (e*^'")+ . On the other 
hand, the product of xSO and phase does not show such behavior. Fig. 4a and 
4b show the real and the imaginary parts of the xSO as a function of (3. We see 
that for small /?(/? = 2.0 and 3.0, for example), the real part has some positive 

15) 

value within error, but becomes unfortunately too noisy to extrapolate to the 
larger j3 region (see Table II) . We should note that the choice of the wave function 
in Fig.l explicitly breaks the parity. This could be the reason why the xSO may be 
non- vanishing for small (3. It is then important to study how the operator behaves 
in the large (3 limit. Since the spontaneous symmetry breaking never happens in a 

24) 

finite volume in the strict sense , the xSO goes to zero as (3 increases. As in the 
case of an introduction of an external source to break explicitly a symmetry, we 
must study finite size effects how the operator becomes vanishing. This is under 
investigation. 
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We have got better results as to the Green's function and the spin-spin cor- 
relations. Fig.5 compares them with the Lanczos method. According to Parola et 
^^15,16) ground states are found to be threefold degenerate; one is the d-wave 
state of momentum (0,0) and the others are a pair of (0,7r) and (7r,0). Both the 
states show slightly different behavior for the density distribution. The Monte 
Carlo results are in good agreement with that for (0,0). Spin-spin correlations are 
shown in Fig.6. The data used for drawing the figures are listed in Table I. 

4-4 Hole doped cases Mf = 5 

The case Mf — 5 makes a "closed shell" in the single free particle distribution 
in the momentum space. The sign problem is very much mild compared to the 
Mf = 7 case. The average of the real part of the phase is shown in Fig. 7 and Table 
II. One sees an exponential fall-off behavior for a wide range of (3 with a decay 
rate being quite small compared to the Mf = 7 case; AE = 0.0018(3). We see 
then that AE depends much on the fiUingness. The unweighted Green's function 
as a function of /? are very stable and gets very small errors. These facts lead us to 
expect that the weighted average is approximately given by the unweighted one. 
In other words, a correction to the unweighted average is quite small. In Fig.8, we 
compare them with the ones obtained by the Lanczos method. We see a very good 
agreement. The spin-spin correlations ( Fig. 9.) are also in very good agreement 
with the exact numerical calculations. The numerical figures are listed in Table 
Ifor the Green's function and the spin-spin correlations. As for the behavior of x, 
both the real and imaginary parts are consistent with zero within error. 
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§5. Summary and Discussion 

We performed Monte Carfo calculations of the two dimensional Hubbard 
model by use of complex wave functions as a trial state. As to the single par- 
ticle density operators, or the Green's functions and spin-spin correlations, we 
have seen good agreements with the Lanczos' results for half-filled and doped 
cases. Encouraged by this fact, we tried to measure the x^O. For half filled case, 
without numerical calculations, we have seen that it is vanishing based upon the 
particle-hole transformation. Upon doping, the phase problem is a main obstacle. 
For Mf = 7, the problem is rather serious, and we are not able to extrapolate to 
large enough values of (3, though the real part of the xSO looks nonvanishing at 
small p. In order to clarify it, we need to study finite size effect to see how this 
operator reaches zero, since no spontaneous symmetry breaking occurs in finite 
volume in the strict sense. This issue is under investigation. For Mf = 5, the 
phase problem is very mild. In this case, the real part of the x^O is consistent 
with zero independent of values of (3. As to the average of plaquettes, which is a 
product of four links (xi2X23X34X4i)) we will report in the forthcoming paper. 
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TABLE CAPTIONS 

Table I The average values and the statistical errors of the single particle density 
operators (spin averaged symmetrical Green's function) Gi and the spin-spin 
correlations Si for various fillings. Suffix i stands for the locations of the sites 
along the path as depicted in the inset of Fig. 2a. The values for Mf = 8, M/ = 
7 and M/ = 5 are shown at /? = 5.0, (3 = 4.0 and (3 = 20.0, respectively. Note 
that Go = 1.0 for Mf = 8. 
Table II The real part of the average of the phase vs. (3 for Mf — 7 and Mf — 5. The 
real and imaginary parts of the xSO is also listed for (3 = 1. Meas stands for 
the number of measurements {Ik = 1000). 



FIGURE CAPTIONS 

Fig. 1 Illustration of the "Dual choice" of wave functions in the momentum space, 
and {il)m} are shown (see text). 

Fig. 2a Expectation value of single particle density operator, or spin averaged sym- 
metrical Green's function Gi for half-filling and {3 = 5.0. The horizontal axis 
stands for the locations i along the path shown in the inset. The first and 
third sites are odd, and the rests are even spacings. Circle is the results of 

16) 

the Monte Carlo calculations, while cross is the Lanczos' results. Lines are 
drawn for the guide of eyes. 
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Fig. 2b Spin-spin correlations Si indicate clear antiferromagnetic correlations for half 
filling and /? = 5.0. The horizontal axis stands for the same location i as in 
Fig. 2a. Monte Carlo calculations ( circle) are in good agreement with those 
of Lanczos' method (cross). The errors lie within the symbols. 

Fig. 3 Average of the real part of the phase of the Monte Carlo weights vs. (3 for 
Mf = 7. 

Fig. 4a Real part of the xSO vs. /3 for Mf = 7. 
Fig. 4b Imaginary part of the x^O vs. (3 for Mf — 7. 

Fig. 5 Single particle density operator Gi for Mf = 7 and /?=4.0. We compare the 
Monte Carlo results (circle) with the Lanczos' ones. According to Parola 
et. al., the degenerate ground states (see text) show different behavior. 
Monte Carlo results agree well with the behavior of the (0,0) state (cross). 
The symbol (square) stands for the behavior of (0,7r) and (7r,0) state. 

Fig. 6 Spin-spin correlations Si for Mf = 7 and /5=4.0. Monte Carlo results (circle) 

16) 

agree with the Lanczos' results ( cross and square). 
Fig. 7 The real part of the phase of the Monte Carlo weights vs. /? for Mf = 5. An 

exponential fall-off behavior is seen in the region where /3 > 10.0. 
Fig. 8 Single particle density operator Gi for Mf = 5 and /9=20.0. The Monte Carlo 

results (circle) are in very good agreement with those of Lanczos' method 

(cross). The errors lie within the symbols. 
Fig. 9 Spin-spin correlations Si (circle) for Mf = 5 and /?=20.0. Comparison with 

16l 

the Lanczos' results (cross). The errors lie within the symbols. 
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APPENDIX 

In this appendix, we show that {PI123 — PI1Z2) appearing in (3.23) is pure 
real. x^O is expressed by two through six point functions. Four point functions 

(Al) 

(1223). ^ {clc2A.cs.) = eg) - G^G^ + Gg^Gg) 
where G^j^ = (cj^Ci.) for spin a. {PI123 - Ph32) in (3.23) 

PZl23 - PI132 = 5Z {(cLc2a4^,C3a'4^„Cl.") - (c}^C3a4a'C2a'4cT"Cla")} 
a,a',a" 

iA.2) 

is written as 

^^^123 - ^^^132 = (112233)t - (113322)| + (1223)|(31)| - (1332)t(21)| 

+ (1231)t(23)i - (1321)^(32)4 + (2331)t(12); - (13)t(3221)4 
+ (I2)t(2331)i - (13)^(3221)^ + (23)^(1231); - (32)1(1321)^ 
+ (31)t(i223)i - (21)^(1332)^ + (il2233)| - (113322)^ 

= — /22/3lfi'l3 + /22/l3fi'31 + 5'22/3l5'l3 — 5'22/l35'31 

+ (terms without f22 or g22) 



(^.3) 



(']) ( I ^ 

where we used the notation = G\ - and gij = G)j . The sum rule (3.21) and 

(A. 3) lead to 

PI123 - PI132 

= {9l39l3 - 93l93l) - (922 + 922){9l39l3 - 93l93l) 

{AA) 

+ (terms without f22 or g22) 

= pure real + (terms without f22 or g22) 
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The terms without /22 or g22 in (A. 4) can also be rewritten by the sum rule (3.21) 
in the form as sum of terms a-\-a* and where a and j3 are complex numbers. 
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